Astronomy & Astrophysics manuscript no. gordon'aa 


©ESO 2012 


April 17, 2012 





Detection noise bias and variance in tlie power spectrum and 
bispectrum in optical interferometry 

J. A. Gordon and D. F. Buscher 

Astrophysics Group, Cavendish Laboratory, University of Cambridge, J.J. Thomson Avenue, Cambridge CB3 OHE, UK e-mail: 
j . gordondmrao . cam . ac . uk 

Received September xx, xxxx; accepted March xx, xxxx 

ABSTRACT 

Context. Long-baseline optical interferometry uses the power spectrum and bispectrum constructs as fundamental observables. Noise 
arising in the detection of the fringe pattern results in both variance and biases in the power spectrum and bispectrum. Previous work 
on correcting the biases and estimating the variances for these quantities typically includes restrictive assumptions about the sampling 
of the interferogram and/or about the relative importance of Poisson and Gaussian noise sources. Until now it has been difficult to 
accurately compensate for systematic biases in data which violates these assumptions. 

Aims. We seek a formalism to allow the construction of bias-free estimators of the bispectrum and power spectrum, and to estimate 
their variances, under less restrictive conditions which include both unevenly-sampled data and measurements affected by a combi- 
nation of noise sources with Poisson and Gaussian statistics. 

Methods. We used a method based on the moments of the noise distributions to derive formulae for the biases introduced to the power 
spectrum and bispectrum when the complex fringe amplitude is derived from an arbitrary linear combination of a set of discrete 
interferogram measurements. 

Results. We have derived formulae for bias-free estimators of the power spectrum and bispectrum which can be used with any linear 
estimator of the fringe complex amplitude. We have demonstrated the importance of bias-free estimators for the case of the detection 
of faint companions (for example exoplanets) using closure phase nulling. We have derived formulae for the variance of the power 
spectrum and have shown how the variance of the bispectrum can be calculated. 

Key words. Instrumentation: interferometers - Techniques: interferometric - Methods: analytical 



■ 1. Introduction 

] Astronomical interferometers combine the light from multiple 

■ telescopes to form interference fringes. The amplitude and phase 

■ of the fringes corresponding to interference between a given pair 
, of telescopes can be combined into a single complex number, the 

■ fringe complex amplitude, which in an ideal interferometer is 
" proportional to a single Fourier component of the object bright- 
. ness distribution. 

" The complex amplitude forms the fundamental observable in 
. phase-stable interferometric experiments, for example at radio 

■ wavelengths, but at optical wavelengths the phase of the fringes 
, is unstable on timescales of order milliseconds, because of op- 
' tical pathlength perturbations introduced by the Earth's atmo- 
sphere. As a result, observers need to combine the information 

. from multiple short-exposure interferograms to recover the de- 
sired science data. Since the phase of the fringes is randomly 
fluctuating, simply averaging the complex amplitude over mul- 
tiple interferograms would lead to a vanishingly small result. 
Instead two observables, the power spectrum and bispectrum, 
are calculated from each interferogram and these can be av- 
eraged over frames as both are immune (to first order) to the 
atmospheric and instrumental phase perturbations. The power 
spectrum serves to capture information about the modulus of 
the complex amplitude, while the phase of the bispectrum (also 
known as the "closure phase") captures information about the 
phase of the complex amplitude on closed triplets of baselines. 

At the low light levels typically encountered in astronomical 
interferometry, the effects of the noise arising in the detection 



process become important. This noise can usually be modelled 
as a mixture of Poisson noise, for example photon noise, and 
signal-independent Gaussian noise, for example readout noise in 
the detector electronics. We use the term "detection noise" to re- 
fer to the combination of all the noise processes in the detection 
process including the photon noise. 

Detection noise produces random measurement errors which 
can be reduced by averaging over many interferograms, but in 
the case of the power spectrum and bispectrum it also produces 
systematic biases in the averaged estimates. There has been con- 
siderable study in the literature on the se biase s (see f o r example 
Goodman & Belsher (1976ab); Sibi lle et al.1 d 19791) : IWirnitzed 
(.1985il : iBeletig (il989D : IPerrinI (12003b ). Some authors have de- 
rived formulae for the power spectrum bias which allow for 
both Poisson and Gaussian noise sources, but existing bispec- 
trum bias estimates have been derived under the assumption of 
Poisson noise only. As a result, it is difficult to derive an un- 
biased estimate of the closure phase from observations where 
both Poisson and Gaussian noise are significant. Such conditions 
are increasingly present in interferometric observations made at 
near-infrared wavelengths so the absence of such an estimator 
is becoming a critical limitation to the acquisition of accurate 
interferometric data as we will demonstrate in Section[8] 

Much of the existing power spectrum and bispectrum anal- 
ysis relies on the direct Fourier transform of the interferogram 
to recover the fringe complex amplitudes. Even under noise-free 
conditions the use of a Fourier-based estimator will return inac- 
curate estimates of the fringe complex amplitude if the interfer- 
ogram is, for example, non-uniformly sampled. To account for 
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this. iTatulli & LeBouquirJ (l2006t) introduced a pixel-to-visibility 
transformation matrix (P2VM) for deriving fringe amplitude es- 
timates in circumstances where the Fourier-based estimators are 
inaccurate, which includes the Fourier estimators as a special 
case. 

There has been considerable study of the variance of t he 
power spectrum and bispectrum (Goodman & Belsher 1976aS 
Daintv & Greenawav 1979; Sibille et al. 1979; Wirnitzer 1981 
Roddier 1987; Avers et alj 1 19881) . but only the work of 
Tatulh & LeBouquin (.2006) applies the more general P2 VM ap- 
proach, and in this case the variance is calculated under assump- 
tions valid only for purely Gaussian noise processes. Knowing 
the variance, and hence the signal-to-noise ratio is important for 
predicting the integration times and detection limits prior to an 
observation run. 

This paper builds on th e introduction bv iGordon et 

(1201 0) and previous work bv iThorsteinss on & Bus cherl (l2004t) 
to present a general formalism for computing biases and vari- 
ances in the presence of detection noise. Our formalism has two 
advantages over current treatments: (a) it can be applied when 
the fringe complex amplitudes are derived using any linear es- 
timator, including the P2VM and Fourier transform as special 
cases and (b) it yields accurate answers for any additive mixture 
of Poisson noise and Gaussian noise. We derive formulae for the 
biases affecting both the power spectrum and bispectrum under 
these conditions and show how our results allow the construction 
of bias-free estimators. We present simulations to show that us- 
ing our new bispectrum estimator can significantly enhance the 
accuracy of measuring the closure phase under realistic observa- 
tion conditions. We also use our formalism to derive expressions 
for the variance of the power spectrum, and explain the method 
by which the reader can calculate the bispectrum variance if re- 
quired. 

2. The Interferometric Equations 

2.1. The Interferogram 

We will start from a specific idealised example and use this to 
introduce our more general model which will be applicable to a 
wide range of realistic scenarios. In our example, the light from 
Mel telescopes is combined to produce a 1 -dimensional interfer- 
ence pattern, which consists of the superposition of a set of sinu- 
soidal fringes at different fringe frequencies /'-' corresponding to 
the interference between telescopes / and j. The interference pat- 
tern is sampled to produce a set of A^pix intensity measurements 
{ip}. We shall call this discrete set of intensity measurements the 
"interferogram". In our example the interferogram is given by: 



^00 1 ^Id 
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H- noise. 



(1) 



where denotes the complex amplitude of the fringes corre- 
sponding to interference between light from telescopes / and j, 
C"" denotes an offset due to the light from all the telescopes (we 
shall call this the "zero-spacing flux"), and is the coordinate 
of pixel p. 

Recognising that Eq. ([TJ represents a sum over a set of 
Fourier components, we can derive an estimate c'' of the com- 
plex amplitudes from the sampled intensity data using a discrete 
Fourier transform (DFT): 



where we note that this equation can also be used to estimate 
the zero-spacing flux by setting = 0. In noise-free condi- 
tions, the estimator c'-' will be equal to the true amplitudes C,y 
provided that the samples are evenly spaced in the coordinate a 
and the frequencies are integer multiples of the fundamental 
frequency given by /fundamental = l/[A^pix(a;2 - a\)]- We wifl call 
this set of conditions the "DFT conditions". 

We can rewrite the above equations in the form: 



and: 



/ = WC + noise, 



Hi, 



(3) 



(4) 



respectively where i — [i\,i2, . . .] is a vector containing the pixel 
intensities, C = [C°°, RefC'^}, Im{C'2|, Re{C'^), Im{C'^), . . .] 
is a vector describing the true complex amplitudes, c - 
[c"", ^{c'^}, Im{c'^), Re{c'^), Im{c'^), . . .] is a vector describing 
the estimated complex amplitudes, W is a matrix known as 
the "visibility-to-pixel matrix" (V2PM) and H is the pixel-to- 
visibility matrix (P2VM) mentioned before. When the DFT con- 
ditions hold, the matrices W and H will consist of sine and co- 
sine terms. 

If we now assume that the DFT conditions do not hold, for 
example if the samples are not evenly spaced or there are "dead" 
pixels, then it is often still possible to set up a linear set of equa- 
tions in the form of Eq. ^ by modifying the matrix W. For the 
case of any "working" beam combiner, we can always find the 
corresponding matrix H such that: 



HW = /, 



(5) 



where / is the identity matrix, ie. Eq. (|5]l is the condition that H 
is a pseudo-inverse of W. It is then clear that Eq. (|4]i can be used 
to derive a value for c such that c - C under noise-free condi- 
tions. The pseudo-inverse generalises the Fourier inverse to the 
cases in which the DFT conditions may not hold; in general the 
pseudo-inverse matrix will not consist of simple Fourier terms. 

Equation (HJi serves to define the basis of our more general 
model for an interferometric measurement. We assume only a 
discrete set of intensity measurements [ip] related to a set of 
complex fringe amplitudes {C). We do not require a simple 
Fourier relationship between the amplitudes and the pixel inten- 
sities; we only require that there exists a set of estimators for 
[Cij] which are linear in [ip]. In other words we require only that 
there exists at least one set of complex weights {//p } such that 
the complex amplitude estimators: 



(6) 



(2) 



have the property that c'-' - C'^ in the limit where detection noise 
can be neglected. This generalised model can be applied to a 
wide variety of beam combination scenarios: the interferograms 
can be sampled temporally or spatially or a combination of the 
two, and multiple telescopes can be combined "all-in-one" into 
a single multiplexed fringe pattern or pairwise onto separate de- 
tectors. 

Note that pseudo-inverse techniques are only presented as 
one possible route to deriving |//j/|; under circumstances where 
the forward matrix W is not known or not fixed (for example if 
the fringe pattern has not been spatially filtered) alternative tech- 
niques to derive {//]/), perhaps based on knowing the statistical 
properties of W, may be needed. Nevertheless, the results pre- 
sented here will be valid independent of the technique used to 
derive values for H'p. 
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2.2. Interferometric observables 



3.3. Other Sources of Noise 



Under phase-unstable conditions, it is common to choose the 
power spectrum, S and bispectrum, B'^'^ as the interferometric 
observables: 



and: 



S'J ^ C'J 



(7) 



(8) 



Typically the values of C'-' will change randomly from interfer- 
ogram to interferogram due to atmospheric disturbances, but the 
values of 5'''' and B'-'* will be more stable. However, unless oth- 
erwise specified, we do not assume any particular distribution 
for the values of Cij; for example the results will apply equally 
well under incoherent conditions where the average value of Cij 
is zero or conditions where a fringe tracker is being used and the 
value of Cij is stable from frame to frame. Our aim is to deter- 
mine suitable estimators, Sq, Bq*, for the power spectrum and 
bispectrum respectively that do not suffer from bias in the pres- 
ence of noise. For this analysis these observables are calculated 
on an interferogram-by-interferogram basis and then averaged. 
Derived quantities, such as the closure phase, are calculated us- 
ing the averaged observables, rather than on an interferogram- 
by-interferogram basis themselves. 



3. Noise Model 

The signal recorded by a real detector system will always have 
some noise component, which becomes especially important at 
low light levels. We model two distinct sources of noise, photon 
shot noise and detector read noise. 



3.1. Photon Noise 

Photon noise is always present in a real interferogram as a result 
of the quantum nature of the detection process which converts 
the incoming photons to an electrical signal. For a given pixel 
p, the probability of A^^ photoevents occurring is given by the 
Poisson distribution: 



a: 



I\p. 



(9) 



where Ap is the ideal noise-free intensity, also called the "classi- 
cal intensity" or "high-light-level intensity". Note that this equa- 
tion implicitly defines the normalisation of Ap in units of photo- 
events, and we require that ip = Ap under noise-free conditions, 
thereby defining the normalisation of ip. 

3.2. Read Noise 

Read noise is present in all optical and IR detectors and is of- 
ten significant at low light levels. Assuming the read noise is 
independent of the number of photoevents, the probability of a 
given measurement error = ip - Np is given by the Gaussian 
distribution: 



^Gaussian ^^pl^/j? ^ 



1 



exp 



(4 



(10) 



where cr^ is the variance of the noise on pixel p; note that the 
variance is not assumed to be the same on all pixels. 



There are additional noise sources that follow the same statistics 
as photon or read noise. Detector dark current. Dp, for example is 
also described by the Poisson distribution and could be included 
in our formalism by replacing Ap with Ap +Dp in Eq. (|9]l. For the 
sake of readability dark current will not be explicitly included in 
this paper 

3.4. Combined Noise Model 

If the photon noise and read noise are statistically independent 
then the probability of obtaining the noisy data, ip, from the un- 
derlying ideal data, Ap, is given by: 

Ppoisson{Np\Ap) 

^Gaussian 

(^ep\o-p,Np)dNp. 

(11) 

The combined noise distribution is a convolution of the Poisson 
noise and Gaussian noise distributions, so it is convenient 
to derive the moments of our noise model via the moment- 
generating function of the probability distribution. The moment- 
generating function for the combined noise model probabil- 
ity density p(/p|Ap,crp) is proportional to the product of the 



moment-generating functions Mpoisson and Mca 
vidual distributions: 



of the indi- 



M (v) = Mpoisson (V) Mcaussian (v) = eXp 



2 2 



exp [-Ap + Ape"] . 

(12) 

The first four moments of the noise model as required to derive 
the results given in this paper are: 













dp) -- 


--^M{v)\ 

dv iv=o 


= Ap 




(13) 


Hi) = 






+ Ap + 0^, 


(14) 


(i,,) = 


d' 1 




+ 3Ap + Ap + 3ApCr^, 


(15) 


at) - 


d" 1 

— M(v) 

dv* lv=0 




+ 6Ap + 7Ap + Ap + 6ApO-p 





+ 6ApCr^ + 3cr^,, 



(16) 



where (. . .) denotes taking the expectation over the probability 
distribution of the detection noise. For simplicity, only the statis- 
tics of the measured data for fixed {Ap} is considered at first; 
averaging over multiple interferograms will be considered later. 
As such the angle brackets do not denote taking an average over 
the distribution of possible values of {Ap}, due to, for example, 
atmospheric distortions to the interferograms. 

An important assumption of the following analysis is that the 
detection noise on different pixels is statistically independent so 
that for example {ipiipi^ - {ipx^iipi^ for p\ p2. It should be 
noted that this assumption will not be true in all cases of prac- 
tical interest. For example when temporally-sampling an inter- 
ferogram with a near-IR detector, it is common not to reset the 
charge accumulation capacitors at each pixel between successive 
time samples, so that a flux measurement ip is derived from the 
difference between the voltage readings on two successive sam- 
ples; in this case the final reading for one sample represents the 
initial reading for the next sample and so the noise on the two 
samples is correlated. Treating such readout schemes is beyond 
the scope of the current paper. 
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4. Bias-free power spectrum estimator 

We seek an estimator S'^ for the noise-free power spectrum S'-' 
that does not suffer from statistical bias introduced by the noise 
on the interfe rogram. A bias-free estimator was calculated by 
iTatulh & Duv ert (2007) in the presence of photon noise and read 
noise, but we present an alternative derivation here to demon- 
strate our analysis method. For the purposes of evaluating the 
power spectrum, we drop all telescope-pair-dependent super- 
scripts, for example writing Hp as Hp. 

We start by evaluating the expectation value of the biased 
estimator |c|^. Combining Eqs. (|6|l and O we get: 



p2 



pi pi 



(17) 



We split the sum in Eq. ( fTTb for the cases where p\ - p2(= p) 
and for p\ p2 which gives: 

p p\ i= p2 

noting that {ip\ ip2^ = (^ip\ ^ (ip2^ for pi p2. Substitution of the 
moments of the noise model given in Eqs. (flSl i and ( fT4l l leads to: 

(kP) = J] {Al + Ap + crl) \Hp\' + Y,^p,Ap2HpiH;,. 

p pi i= p2 

(19) 

The noise-free power spectrum, S , can be expressed in terms 
of the noise-free interferogram, Ap, as: 



YYApiAp2HpiH*p2' 

pi p2 



(20) 



This can be split into terms equivalent to those in Eq. ( fT9T l: 

^ = Z ^p 1^/ + Z T,^p^^p^"p^"p2' 

p pi + pi 

which can be substituted into Eq. (fT9] l to give: 
(|c|2) = 5+2A„|//„f , 



(21) 



(22) 



where Nf denotes the number of frames over which the average 
is taken. This result depends on the assumption that the detection 
noise is dependent only on the value of the flux A^ in each pixel 
and is statistically independent of the particular shape of the in- 
terferogram and hence of the atmospheric distortions affecting a 
particular interferogram. 

For the case where read noise is negligible (a-p — Oe" and 

H is the DFT (\Hp^ ^_|e^™W]^_^l) Eq. dH reduces to the 
estimator reported bv lGoodman & Bels her ( 1976a) : 

5i=|c|2-A?, (26) 

where N = 2p is the total number of photons in an interfer- 
ogram. For high light levels (negligible photon noise and read 
noise), the \cf term dominates (as this goes as N^), and we re- 
duce our estimator further to the uncorrected estimator (correct 
only under zero noise conditions): 

.|2 



Si = \c\' 



5. Bias-free power spectrum variance 



(27) 



The derivation for the power spectrum variance requires the first 
four moments of the noise model, Eqs. (fTST i to (fTSI l. and stretches 
to multiple pages due to the \c'^'^ term. The derivation follows 
the methodology presented above and in the derivation of the 
bias-free bispectrum estimator given in Appendix |A] but is not 
presented in this paper For readability, the formulation of the 
power spectrum variance, var(5o), is split into terms that dom- 
inate in the case where the noise is dominated by photon noise, 
read noise or where both noise sources are significant: 



var(5o) = 



p 







|2 









'-{Sf 



Photon Noise Terms: 

|2 



= 2 |C|2 l^pf + <^<^ Z ^pKK + Z ^p"p"p 

p p p 

r 1- 

+ Z^'^K^I' +Y,^pHpHpY,^pH;H; 



and written in terms of the noisy interferograms we see: 

|2\ 



Read Noise Terms: 



p 



The bias-free power spectrum estimator, Sq, is then: 

p 



>o 



ip + (t]) \Hp\ 



(23) 



(24) 



It needs to be remembered that is intended to be evalu- 
ated on a per-interferogram, ie. frame-by-frame basis as atmo- 
spheric turbulence changes the shape of the interferogram be- 
tween frames. It is straightforward to show that when this esti- 
mator is averaged over a number of frames that the average is 
itself unbiased. This is because taking the average over frames 
commutes with taking the expectation over the detection statis- 
tics, so that: 



P L p 

Coupled Terms: 



+ 2 2 Yu "^p \"p\^ + Z ^p"p"p Z ""l^i'^'p 

p p p p 

+ Yj ^'KK Z "^p^p^p + 2 Z "^P 1^/ 
p p p 

+ ccY, orlH;H; + cc ^ cj^pHpHp -icYo-\ \hJ( h; 

p p p 

- 2C'Y,^l\HpfHp. 

p 

For comparison with previous work we rewrite Eq. ( l28T l assum- 
ing DFT conditions, that read noise is constant across all pixels 
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10 10 
Number of Source Photons on Detector 



Fig. 1. The relative importance of the photon, read (cr^ = 3e ) 
and coupled noise terms under DFT conditions for the power 
spectrum variance with respect to the number of source photons 
detected. The area highlighted in grey is the region where the 
coupled terms dominate the power spectrum variance. 



(Hp - Npix<T^), and that atmospheric phase fluctuations are 
random and as such the cross-product of independent complex 
quantities average to zero: 

Photon Noise Terms 



var (5dft) = 2N \C\^ +N^ + Y, ApC^'^'^-'f'" ^ Ape-*"''""^" 



Read Noise Terms 



Coupled Terms 



+ 2\C\^Np,,cr^ +2NNpi,cr\ 



(28) 



Separate analysis of t he photon and detector no i se terms, 
as is o ften done (see e.g. 'Perri n & Ridgwavl (l2005h : iBouquinI 
means that the Gaussian-Poisson coupled terms are 
missed. Inspection of Eq. ( |28] | and Fig.[T]show the importance of 
including these coupled terms for a wide range of photon fluxes, 
and therefore of incorporating both the photon noise and read 
noise into a single noise model. In the above example, for ~ 300 
source photons on the detector, the variance on the powerspec- 
trum is underestimated by a factor of over 6 if cross-terms are 
excluded from the variance calculation. 



lengthy and is given in Appendix |A] We present here the bias- 
free bispectrum estimator, Bq*: 



= c''c^''c^' 
- c' 



p 

p 
p 

+ Y,{'2ip + 3cTl)H'jHfH'. 



(30) 



This equation has been verified though computational simula- 
tions to the 0. 1% level. In sections|2]and[8]we present the results 
of these simulations regarding the closure phase, derived from 
the bias-free bispectrum estimator in Eq. (l30l l. 

If we assume DFT conditions so that Hp = e^™pf'' and and 
an all-in-one interference pattern such that + /-'* + f'^' - 0, 
then we have Hfu'^p - Wp where we note Hjl - (^H'J^ and 

H'jHfH^ = 1. 

If we also assume the read noise cr^ = cr is constant across 
all pixels rather than varying on a pixel by pixel basis then: 

^o-V™'"^"=0. (31) 



The fact that the third order moment of the Gaussian distribu- 
tion is zero would suggest that there should be no read noise 
terms under DFT conditions for constant cr, however it must be 
remembered that there are second order moments in the bias cor- 
rection terms which will contribute to the bias under these con- 
ditions. 

If we further assume that read noise is negligible (o-p - 0), 

then Bq* reduces to the estimator presented bv ^ Wirnitzerl (1 1 985h 
accounting for the bias introduced by photon noise alone: 



■2N, 



(32) 



where - Yjp ip is the mean number of photons in the interfer- 
ogram. 

For high light levels (negligible photon noise and read noise), 
the c'^d'^c'^' term dominates (as this goes as N^), and we reduce 
our estimator further to the uncorrected estimator (correct only 
under zero noise conditions): 



6. Bias-free bispectrum estimator 

The bias-free bispectrum estimator is calculated in the same 
manner as the bias-free power spectrum estimator. The defini- 
tion of the bispectrum presented in Eq. ([8]l expands to: 

B'V* = C^d^c'^ = Z Z Z VA,3//X2<3- (29) 

p\ p2 p3 

Splitting the sum in Eq. ( |29] l gives 5 terms rather than 2 as for 
the power spectrum estimator. As such the derivation is quite 



fif = c'V^c'^'. (33) 

The bispectrum variance is not presented in this paper as the 
\c'^c^''c^'\^ term leads to a long and tedious derivation. The zeal- 
ous reader is invited to calculate the expression for the bispec- 
trum variance if required using the method presented for calcu- 
lating the bispectrum in AppendixlAl 

The remainder of this paper will be concerned with demon- 
strating the magnitude of the bias introduced when the wrong 
estimator is used in the presence of photon and read noise at low 
light levels. 
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-B;-" ,-B;' Estimators 
k;-* Estimator 
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Object Closure Phase [degrees] 




— Jj;'' Estimator 

- - J3* Estimator 
B'-* Estimator 



-60 60 

Object Closure Phase [degrees] 



120 




-180 -120 -60 60 120 180 -180 -120 -60 60 120 180 



Object Closure Phase [degrees] Object Closure Phase [degrees] 

Fig. 3. The bias in closure phase for By * (soHd Une), B j'* (dashed Une), and B'^'' (dotted Une) estimators under a number of measure- 
ment conditions. Plots a and b show DFT conditions, whilst plots c and d show non-DFT conditions (described in the text). Read 
noise of o-p - 3e" is present in plots b and d. The thin and think lines show 120 and 300 photons reaching the detector respectively. 
The shaded regions in plots c and d show the range of possible biases resulting from different possible combinations of phases on 
the three baselines contributing to the source closure phase. 



7. Closure phase bias 

To illustrate the magnitude of the systematic bias in a realistic 
observation we present the closure phase data from a set of sim- 
ple simulations. We created ideal interferograms resulting from 
a three telescope interferometer observing a target with a known 
and controllable closure phase. Photon noise and read noise were 
added, and the bispectrum (and hence closure phase) extracted 
using each of the three estimators presented above. The interfer- 
ograms were created under either DFT or non-DFT conditions. 
The fringe frequency on the detector for each baseline was set 
at a ratio {1:2:3) for the DFT case, and 1.1 x {1 : 2 : 3) for 
the non-DFT conditions. To simulate non-DFT conditions we 
rewrite Eq. ([T]) as: 



Mel 

ip = bpC'^ + 2 ^ + noise, (34) 

•<j 



where the beam profile, bp is either flat: bp - I /Np\^ for DFT 
conditions (as in Eq. ([1])), or a cropped Gaussian of the form: 



exp 




V- 




i:?"exp 


(p'- 

2f 


2 



(35) 



where ^ = 22 pixels and fi is centered on the 64 pixel de- 
tector From Eq. ( l34l l we can construct a V2PM, W, which is 
correct under DFT and non-DFT conditions. The form of W is 
shown in Fig. |2] for the non-DFT conditions described. We use 
the Singular Value Decomposition method to find the optimum 
(in terms of least squares) pseudo-inverse of W, which gives us 
H. An instrumental visibility loss was introduced such that the 
maximum visibility on a baseline is |y,| - 1/3. The results pre- 
sented in Fig.[3]are calculated using the ideal interferograms and 
therefore show only the closure phase bias and suffer no variance 
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u ■ • 

Q - : 

4 8 12 16 20 24 28 32 36 40 44 48 52 56 60 
Pixel Number 

Fig. 2. The transform matrix, W under non-DFT conditions. The 
DC term shows the beam profile adopted, and the cosine (soUd 
hnes) and sine (dotted Hnes) components of the fringe pattern for 
each basehne show the deviation from integer number of fringe 
periods across the detector (with frequencies in the ratio 1 . 1 x { 1 ; 
2 : 3)). This is similar (but with the addi tion of a DC term) to 
the vis ibility to pixel matrix introduced bv lLe Bouquin & TatulUl 
(l2Q06h . 

contribution. These are given by: 

Bf = B-j' + A.e^''-''/'^ 2 (a + cr^) e^-""/^' 
p p 
+ ^ A„e2--"/^' ^ (a + cr^) e^-"^^'^' 

p p 
+ Ape^''*""/" ^ (a + e^'''""f"' + n, (36) 
p p 

and: 

gijk ^ gUk _ g ij _ g jk _ski_j^_ 2Npi^O-^ , (37) 

whilst the bias free estimator is simply: 

B'f = B'^K (38) 

We note the order of the indices in the DFT terms in Eqs. ( |36] | 
and (l37T i and that the ideal source bispectrum and power spec- 
trum are independent of the measurement system (weather or not 
if meets the DFT conditions), whereas the bias terms implicitly 
assume DFT conditions for the Bj'* and Bj* estimators. Fig. [3] 
shows a simulation of the bias present on the closure phase for 
the three estimators for simulations of 120 and 300 photons per 
interferogram reaching the detector. 

7.1. Zero read noise 

In the situation where photon noise dominates the noise model 
((Tp - Oe") and the instrument is setup such that the interfero- 
grams are formed under DFT conditions, Bq* and Bj'* both re- 
turn the correct closure phase, however Bj* is affected by a bias 
term resulting from photon noise. At these light levels, the bias 
in Bj* becomes significant for the majority of source closure 
phases. 



7.2. Non-zero read noise 

At low light levels the read noise on the detector becomes im- 
portant and comparable to the photon noise. Fig.[3j) shows sim- 
ulations under DFT conditions, as in Fig.|3^, but with o-p - 3e" 
read noise included. This level of read noise is quite low com- 
pared to many of the optical detectors in use today. We see that 
failing to account for read noise introduces a significant bias in 
closure phase estimates for both B^* and B'^''. It is clear that the 
magnitude of the bias term is strongly dependent on the source 
brightness, and that Bj * is preferable to B'^'' under certain com- 
binations of noise and photon flux. 

7.3. Non-DFT Conditions 

Deviation from DFT conditions changes the way statistical bias 
affects the bispectrum. Inspection of Eq. ( l30l l shows that the sta- 
tistical bias can affect both the real and imaginary parts of the 
bispectrum when H is not the DFT. If the measurement system 
conforms to DFT conditions, and the read noise is negligible or 
constant across all pixels, the bias only affects the real part of 
the bispectrum estimators as shown when applying Eq. ( |3TI ) to 
Eq. ( l30b . It is also apparent that the magnitude of the bias is 
no longer dependent only on the value of the closure phase, but 
on the value of the individual phase of each baseline. The non- 
uniform beam profile causes an overlap of the fringe compo- 
nents in frequency space, and the non-integral fringe frequency 
on the detector invaUdates the assumption that 7/^*//*' - and 

H'jHfHp = 1 which led to the formulation of Bj''^. Fig. [3]c is 
provided as an example of the impact of non-DFT conditions on 
the closure phase bias for a photon noise dominated instrument 
(cTp - 0). Note that the non-DFT conditions change the bias 

terms, meaning Bj'* is now biased when it was not under DFT 
conditions. Fig. |3]d shows the effect of non-DFT conditions on 
a system with photon and read noise and again there is an al- 
teration of the bias terms which Bq corrects but the other esti- 
mators do not. The grey regions on the non-DFT plots show the 
range of biases associated with each estimator depending on the 
individual phases of the triplet producing a given closure phase. 
The lines within the shaded regions represent the expectation of 
the closure phase bias over all possible individual phases. 

We see that the bias on the closure phase does not cancel 
when accumulating frames; specifically, that there is no advan- 
tage to observing in phase-unstable conditions over phase-stable 
conditions as a means of reducing bias. Therefore, bias-wise, 
one should always favor fringe-tracker assisted observations 
whereby the longer discrete integration times result in higher 
SNR, even in the photon-rich regime where fringe tracking is 
not strictly necessary. 

7.4. Closure pliase bias and statistical uncertainty 

The bias is only significant at low light levels, and tends to zero 
as the light level increases. At low light levels the variance is also 
high, and we must determine the relative importance of these 
two processes. We present Fig. |4] which shows the relative im- 
portance of bias, Ot'i^s, and closure phase standard deviation, g^, 
using different estimators. 

We show the number of frames, Nf, such that: 
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Table 1. Simulation Parameters 




10^ 10^ 10' 10= 

Number of Source Photons on Detector 



Fig. 4. The number of integrated frames required to reduce the 
statistical uncertainty to the level of the bias for the B j'* (dashed) 
and B'J'' (dotted) estimators. 



Parameter 


Value 


Primary Raduis 




Primary K-band Magnitude 


5.0 Mag 


Companion Separation (s) 


(10, 100, 1000) i?* 


Companion flux ratio (r) 


1% 



31 



Telescope Diameter 1.8 m 

Number of Apertures 3 

Instrumental Baselines {1:2 

Optical Throughput 1 % 

Instrumental Visibility 1/3 

Center Wavelength 2.21 lyum 

Bandwidth 48 nm 

Exposure Time 25 ms 

Number of Frames 1 0^ 1 0^ 1 0^ 

Total Integration Time 25s, 4ml0s, 41m40s 

Read Noise cr = Se" 

Quantum Efficiency 80% 

Science Detector Pixels 64 



Notes. The top and bottom panels show the source and instrument pa- 
rameters respectively. 



which is the number of interferograms required to reduce the 
magnitude of the statistical uncertainty to that of the bias. The 
closure phase standard deviation per frame, ^-o, is given by: 



1 



Im 



Bf X 



\B'M 



(40) 



where B'j'^ is the bispectrum calculated per frame using the cho- 
sen estimator and Im ^B'j'' x jfr^ j gives the projection of B'j'' on 

a unit vector perpendicular to B''*, and. This method is used be- 
cause it is insensitive to wrapping errors in the complex plane 
when the variance is large. In this example, the magnitude of the 
bias and variance are measured when the source closure phase is 
at 90°. The values presented in Fig.|4]were derived using numer- 
ical simulations to the 0.1% level, and the bias contribution was 
verified analytically using Eqs. ( |36] | and ( |37] |. 

It can be seen in Fig. |4]that the bias is important for real- 
istic numbers of frames (e.g. 10^ frames or 25s integration for 
25ms exposures) compared to the uncertainty for a wide range 
of typical light levels. 

8. Simulation: Detection of faint companions 

Closure phase nulling is an observ ational techniques that allows 
the detection of faint companions (Chelli et al."2009). The tech- 
nique has been demonstrated by Duvert et al. (2010), detecting 
a 5 magnitude fainter companion using AMBER/VLTI around 
HD 59717. 

Closure phase nulling uses the fact that the closure phase 
jumps from to 180 degrees as the visibility function of the 
primary in the system passes through a null. The presence of 
a companion becomes evident as a deviation from an instanta- 
neous phase jump. The profile of the phase jump gives informa- 
tion about the projected separation and flux ratio of the system 
components. We show that statistical bias in the closure phase 
leads directly to incorrect e stimation of the system properties. 

We follow IChelli etal] (2009) and simulate a binary system 
where the primary is described by a uniform disk of radius, 



with a point source companion at a separation, s, with flux ratio, 
r. For simplicity we simulate a three telescope interferometer 
such that the position angle of the binary system is parallel to the 
frequency axis of the interferometer. The interferometer samples 
three frequencies {mi2, "23, wn) such that = ui2 + M23- The 
instrumental parameters are given in Table[T]and we assume the 
instrument is functioning under DFT conditions. 

The complex visibility, OHu), at a given frequency, u, is de- 
scribed by: 

V.(.)..e-^-^ 
1 + r 

which is the visibility function, y*(M), of the primary: 

J\{2nuR-i,) 



V*(m) = 2- 



(42) 



plus the contribution of the secondary with a phase modulation 
related to the projected separation, s, and flux ratio, r. J\ is the 
first order Bessel function and R^. is the radius of the primary. 

The ideal bispectrum is extracted from Eq. (HTt using Eq. ^ 
and the closure phase is taken as the argument of the averaged 
bispectrum. The closure phase signal around the fi rst null in the 
visibility profile is shown in Fig. 2 of .Chelli et alJ (12009) and re- 
produced using our simulation model in Fig.|5^ in this paper for 
a number of companion separations (s - 10/?*, 100/?*, 1000/?*). 
The closure phase is wrapped modulo 2;7r and the frequencies 
{mi2, M23, M13) are set to the ratio {1 : 2 : 3} on the detector and 
scaled by the radius of the primary, /?*. For simplicity, the base- 
line to fringe-frequency mapping is chosen such that it preserves 
the ratios of the two. The specifics of the mapping is not expected 
to have a significant impact on the bias. The closure phase of the 
primary alone (single source) is shown for comparison. 

We present in Fig.|5]3 a simulation of the closure phase from 
a 100/?* separation binary with flux ratio 0.01. Bq* returns the 
correct closure phase signal in agreement with the ideal case 
presented bv IChelH et al] (|2009|) . fif and Bf are strongly af- 
fected by bias terms, meaning the closure phase is significantly 
different. Depending which regions of the closure phase trace 
are sampled by the interferometer, B'^^ and B'f^ would return an 
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Fig. 5. Demonstration of the effect of bias on a closure phase measurement. Plot a: The closure phase trace resulting from a binary 
pair consisting of a resolved primary of radius and an unresolved secondary with flux ratio r - 0.01 for separations (s - 
10/?*, IQQRi,, 1000/?*). Plot b: The closure phase trace for the 100/?* system recovered by the three estimators. Bq*^ correctly 
returns the ideal closure phase, whilst Bj^* and Bj * do not. Both plots show the single source closure phase signal (thin solid line) 
for reference. Plots c and d show the closure phase calculated using 5^* and * respectively as in plot b along with the statistical 
uncertainty on the measurement after 10^, 10"* and 10^ frames. 



incorrect value for the primary radius and companion flux ratio. 
Figs.jSj; and|5}l show the closure phase calculated using Bf and 
/?2 * respectively and the statistical uncertainty on the measure- 
ment after 10^, 10^ and 10^ frames calculated using Eq. ( l40b . We 
see that even for very short total integration times (lO-' frames is 
under a minute using typical discrete integration times) the bias 
still produces a significant contribution. 

9. Conclusions 

Interferometry offers exceptional angular resolution at the ex- 
pense of sensitivity, as the collecting area is, by necessity, 
smaller than that of a single dish with the same resolution. In 
ground-based optical interferometry these sensitivity limitations 
are accentuated by the need for short observations to reduce the 
phase perturbations introduced by the atmosphere. The result is 
that noise levels on individual observation frames are often high. 



exacerbated by the need for fast read-out of the detector. Without 
understanding the statistical properties of the noise, bias terms 
are introduced when extracting information from the interfero- 
gram using the non-linear power spectrum and bispectrum con- 
structs. We have shown in this paper that: 

1 . Neglecting the cross terms that arise in the presence of both 
read noise and photon noise leads to bias on the extraction 
of interferometric observables resulting from the statistics of 
the noise. 

2. The statistical bias on the observables results in incorrect 
phisical parameters being extracted from observation data, 
for example the flux ratio of binary pairs or the diameter of 
a target as shown for the case of closure phase nulling. 

To remedy this we have presented bias-free estimators for 
both the power spectrum and bispectrum (and hence the closure 
phase) which can be used in the presence of any combination 
of photon and read noise under any fringe encoding scheme in 



9 



J. A. Gordon and D. F. Buscher: Bias and variance in the power spectrum and bispectrum 

which the measurement equation is linear. Our bias-free estima- 
tors are presented in a format allowing easy implementation into 
current data reduction pipeUnes. 

Acknowledgements. We thank the anonymous referee for constructive com- 
ments and insiteful suggestions that have undoubtedly improved this paper. 
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Appendix A: Derivation of the bias-free bispectrum Substituting the ideal bispectrum split sum from Eq. (IA.2b gives: 
estimator 



The ideal, noise-free bispectrum is given by: 

p 1 p2 p3 

which can be spHt into: 

p 

p * pi 



P* p2 



P*pi 



ki 



(A.2) 



The equivalent expression for real data embedded with noise 
when average over frames is: 



p 

p* pi 

p ^ p2 

+ Z Z (^^ + "^D \iH^jHfH% . (A.6) 

pi=p3 

By rearranging and substituting in for the ensemble average over 
the real interferogram we find the ideal bispectrum in terms of 
the noisy interferogram: 

B'j'' = {c'^c^'^c'") 

^ p 

p * pi 



(A.7) 



P * pi 



(c'^c^V*') = Z Z Zl ('pi'/'2'/'3) Ki^p2^pi' ^^-^^ Expanding out the double summations using the fact that = 



p 1 p2 p3 

which splits into: 

(c'^c^^c^') = J]{i^,)H^jHfH^' 

- i:h{^p){^p^)"'j"'>p 

p It pi 

^ ZZZ('^'>)('.^)('.3)^1<2<3- (A.4) 

pi # p2 ,t p3 

Substitution of the noise model in Eqs. ( fTSl l to ( fTsT i gives the 
noisy bispe 
ated noise: 

(c'W') = ^(Af, + 3A2+A, + 3A,(r2)//;/<i/*' 

p # pi 

p # p2 

+ Z Z (^^ + + '^p) ^p^^p^fKi 

p* pi 

+ ZZZVAp2Ap3^</ii^^X'3. (A.5) 



22-2 gives: 

a b a=b 



[j] (3/2 - + 3;X - ^o-])H'lHfH\; 



+ 3 + 



pi 



P2 



- Z('^ + ^.)^'<Z'"3<3 ■ 



^3 



(A.8) 



""'f^^i'P^"^™™ °^ *^ ^^^""^ interferogram and associ- Cancelling terms and noting that, Y.p ipH^^ = c'\ we find the 

bias free bispectrum estimator: 



bias-iree 



- d'Y^^P + ^nfH^ 

p 

p 
p 

+ ;^(2/, + 3cr^)/^^<H^'- 



(A.9) 



pi # p2 ^ p3 



